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Abstract 

With the aim of stndying the relevance and properties of critical slowing 
down in Monte Carlo simnlations of lattice quantnm held theories we car¬ 
ried ont a high precision nnmerical study of the discretised two-dimensional 
QpAT-i at A = 10 using an over-heat bath algorithm. We identify crit¬ 

ical slowing down in terms of slowly-evolving topological modes and present 
evidence that other observables couple to these slow modes. This coupling 
is found to reduce however as we increase the physical volume in which we 
simulate. 

Keywords: Critical Slowing Down, CP^“\ Over-heat bath. Topological 
charge 


1. Introduction 

Two-dimensional CP^“^ models are of great interest to lattice QCD prac¬ 
titioners, as they exhibit similar properties to QCD, such as conhnement and 
asymptotic freedom mi- Because of their lower dimensionality, they are 
also much easier and less costly to simulate; for this reason they provide an 
excellent framework in which to study and test the properties of algorithms. 
For example, they are very useful for studying the phenomenon of critical 
slowing down (CSD): an increase in Monte Carlo relaxation time for particu¬ 
lar observables when generating a new statistically independent conhguration 
on which that observable can be measured. 

Within the context of the CP^“^ model, an exponential form of CSD for 
the topological modes has been well documented |3HS]. Other observables. 




such as the magnetic susceptibility, are expected to largely decouple from 
the topological modes |1] and thus the CSD is expected to conform to a 
simple power law. However potential small deviations from this behaviour 
have been observed at high values of the correlation length [6]. Our aim is 
to study CSD in the model to a uniquely high level of precision, with 

the aim of quantifying any deviations from the expected form of scaling for 
very hne lattices. 

Our detailed studies of CSD require the generation of hundreds of millions 
of lattice conhgurations; as such we are able to determine certain observables 
in the CP^“^ model to a high level of statistical precision. For this reason 
it becomes necessary to be alert to numerical errors that algorithms may in¬ 
troduce, which are usually vastly outweighed by statistical errors. We there¬ 
fore also present the hne-tuning necessary to the over-heat bath algorithm 
employed to overcome such numerical errors that may become manifest in 
particularly long Monte Carlo simulations. 

The layout of this paper is as follows: in section we introduce the 
continuum formulation of the CP^“^ model. In section we provide the 
details of the lattice formulation that we use and detail the observables that 
we measure in our simulations. In section we give an overview of the over¬ 
heat bath algorithm that we use in our simulations, as well as detailing the 
adjustments necessary in order to avoid numerical errors. We then go on to 
present the results of our numerical simulations. In section we hrst give 
an overview of the hnite volume effects of our model, which motivates the 
choice of lattice sizes used in our simulations of CSD, which we report in 
section We demonstrate a deviation from the expected Gaussian scaling 
for the integrated autocorrelation time of the magnetic susceptibility, which 
we can clearly attribute to a coupling between the topological modes and 
other observables. In section we present our conclusions. 

2. The CP^-i Model 

The CP^“^ model in 2 dimensions is classihed by the action nizi: 

S = i y d^x (x) ■ D'z (x) j , (1) 

where ^ is an Wcomponent vector of complex scalar helds subject to the 
constraint 

z{x) ■ z {x) = 1. (2) 


2 



The covariant derivative is defined as 


Dfj, — d^ + ( 3 ) 

and Afj^ is an auxiliary U (1) gauge field. 

Using the composite operator [3] 

P (x) = z (x) z{x), (4) 

we can obtain the correlation function 

Gp (x) =Tr {P (x) P (0))„„„ = Tt (P (x) P (0)> - i (5) 

This function will be of importance for deriving other lattice observables, 
which we will discuss in the next section. 

The main quantity of interest for our simulations is the topological 

charge of the field (x), which is defined by DEI- 

q (x) = = Tepp-DpZ (x) ■ D^z (x). (6) 

ZTT ZTT 

This quantity is related to the topological susceptibility through 

Xt = j d‘^x{q{x)q{0)). ( 7 ) 

3. Lattice Formulation 

The lattice formulation of the model that we use follows from 

employing a first-order discretisation of Eq. Q using a periodic square lattice 
of side length L. This leads us to the action [31 [M] 

Sg N(3 ^ ^ T ZnZxi-\-gXn,g, 2^ , (8) 

n,fi 

where we have introduced the gauge links G U (1) connecting adjacent 
lattice sites and N(3 = 1/g. We index the lattice with with n = {nx,ny) and 
the direction of the link with fi = 1,2. While it is possible to integrate out 
the f/(l) gauge field using its equation of motion, it is useful to leave the 
fields in explicitly. This makes the resulting lattice formulation of the action 
linear with respect to each variable, which allows the use of local updating 
algorithms, such as the over heat bath algorithm. 
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3.1. Lattice Observables 

We can derive two important observables from the Fonrier transform of 
the correlation fnnction, dehned as mu 

^ ^ exp (n - m) ■ fcV (9) 

n,m ^ ' 


where n and m index the points on the lattice, and the 2-momentnm of onr 
theory is dehned on a periodic lattice by p = {2ti/L) k. In a hnite volnme, 
the entries of k only take integer valnes in the range [0,A^). Using Eq. @ 
we can extract the magnetic snsceptibility, Xm, from 

Xm = G(0,0), (10) 


and a dehnition of the correlation length, from 


1 ( g(o,o) _ 

4sin^ (tt/L) yc (1,0) j 


( 11 ) 


We will however make use of another dehnition of the correlation length, 
dehned through the wall-wall correlation function, dehned as |3] 


G„ (n^ - m^) = - ^ Gp (n^, n,; m,). 


( 12 ) 


ny^rriy 


On a lattice with periodic boundary conditions, the long distance be¬ 
haviour of the correlation function takes the form 


(^) ^ ^ 





(13) 


We can therefore extract the parameters and by htting the function 
Gw (x) in a region dominated by the ground state contribution. 

We use the geometrical dehnition of the topological charge, given by |7] 

Qn = {In [TrPn+fp+uPn+f,Pn] + In [^TPn+uPn+^l+uPn]} , (14) 

ZTT 

with the constraint /i 7^ z/. The topological susceptibility of the conhguration 
is then obtained from 

( 15 ) 
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3.2. Autocorrelation 

In order to quantify the CSD of observables in our simulation, we measure 
the integrated autocorrelation time HD]. We start by defining the autocorre¬ 
lation function for an observable O as: 

n=l 

where iVc/ is the total number of conhgurations we have in our Monte Carlo 
sample. The dehnition of the integrated autocorrelation time for this observ¬ 
able is 



C{t) 
^C(O) 


E 


(17) 


It is important to take N^f tq in order to obtain an accurate measure¬ 
ment of the integrated autocorrelation time of an observable. For t > tq we 
begin to capture more noise than signal in our summation, and so in practice 
we employ a windowing function A (t) such that \(t) ~ 1 for t < tq and 
A (f) ~ 0 for t^ To [in] . In our simulations we use the common choice of 


Xit) 


1 N < M, 

0 \t\ > M, 


(18) 


where the cutoff M is chosen such that N^f M > ctq. This is referred to 
as the “automatic windowing” algorithm. This method of course introduces 
a systematic error because of the truncation of the sum. We must therefore 
choose a large enough window in order to keep this error small. Our primary 
source of error on the measurement should therefore be statistical, which we 
can approximate using the estimator uni: 


al ~ 


2 (2M + 1) 


N, 


cf 


IQ. 


(19) 


We also verihed that we obtain compatible results using the automatic win¬ 
dowing method of Ref. m! 
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4. Over-heat bath Algorithm 

In order to undertake simulations of the model we have im¬ 

plemented an over-heat bath algorithm na, which is as follows (following 
Ref. 0). We hrst express the action at a particular lattice point n as a 
scalar product of two vectors. The local action relating to the helds Zn can 
be written as 


Sn,z = -(3NRe {Zn ■ Rz,n} , 

(20) 

and for the helds \n,^ as 


Sn,\ = -/3NRe{Xn,^Fx,n,^,} ■ 

(21) 

The explicit forms of the F terms are therefore: 



(22) 

F 

(23) 


These expressions can be written in terms of real vectors, 0 and with 2k 
components (e.g. arranged such that (j) 2 j = and 02 i+i = Inijznj} 

for i = 0, ...,A; — 1). We have k = N and k = 1 for z and A updates 
respectively. The contribution of the vector 0 to the action is then given by 
(dropping the n subscript for readability): 

S^ = -PN(f)- F^ = -PN\F^\cose. (24) 

To update the vector 0 we simply generate a new angle 9new from the prob¬ 
ability distribution 

= (sin6')^''"^exp(0A|F0|cos6'). (25) 

d cos 0 

The condition for the over-heat-bath algorithm is that the new vector (pnew 
is chosen by minimising the scalar product between (pnew and <poid, which is 
satished by taking 

1 _ n -^0 f j. n ^4’ ^ 9new (oa\ 

\F^\ \ sin6'o/d 

cos sin 
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Our code is vectorized as follows: first we define the as the number 
of doubles that will £t into the vector register of the CPU we employ. We 
employ red-black ordering as a means to parallelize our code; however we 
also group adjacent sites of the same colour (i.e. red or black) into blocks of 
size Nd- We then use vector instructions to generate new angles to update 
these Nd sites simultaneously. By employing this method of vectorization we 
must impose the constraint that L must be a multiple of 2Nd. Running on 
a single Intel Xeon E5-2670 processor, for a single thread we obtain a peak 
performance of ~ 5.2 GFLOPS for the over-heat bath update procedure 
alone (~ 25% of peak processor performance) and ~ 5.9 GFLOPS for the 
entire simulation (~ 28.5% of peak processor performance). On a lattice of 
size L = 160 this corresponds to generating and taking measurements for 
188,500 configurations per core hour. 

Additionally we remark that for all our runs we took care to ensure that 
our Monte Carlo time series are thermalised, by verifying that the average of 
each observable is stable with respect to the number of thermalisation steps. 


4.1. Trial Angle Generation 

The details of generating a new angle according to the distribution Fq. 


(25) are given in the appendix of Ref. [3], which we repeat here for conve¬ 
nience. We cannot directly generate angles according to Fq. (25), and so 


we instead generate angles according to a similar distribution and employ an 
accept/reject step to ensure that the angles we generate do £t the required 
distribution. First a trial variable, 6, is generated according to the Lorentzian 
distribution 


pI (0) = 


l + c^9- 9o) 


2 • 


(28) 


Dehning 


c = 


k-1 


PN\F^y 

We dehne the two parameters of this distribution 9o and c as 


9o = arccos ^V^lT^ - , 

c=\//?iv|fiiyi + e. 


(29) 

(30) 

(31) 
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Figure 4.1: Plots of ry — 1 against PN \F^\ for (a) fc = 1 and (6) fc = 2,4,6,10. 
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The trial variable fitting the distribution in Eq. (28) are obtained from a 
uniform distribution in x G [0,1] using 


9 = 9o -\— tan [x arctan c (tt — 9q) + (x — 1) arctan c6 *o]. (32) 


To ensure that the trial variables generated match our desired distribution, 
we accept it with the probability 


Pk (g) 1 + c^ i9-9o)" 
Pk (9o) P 


(33) 


where pk {9) is dehned as 

Pk {9) = (sin exp (/3iV \F^\ cos 9) , (34) 


and the parameter p is chosen such that Pace < 1- 

In order to choose this free parameter p, we initially set p = 1 and 
determine the maximum value that Pace for a range of values of /3N\F^\. 
The measured values of then represent the optimal choice for p. In 

Fig. |4.1| we display our plots of p for various values of k. For the z updating 
{k > 2), it is clear that this choice of p is only somewhat of a concern for 
very low N. While p = 1.01 would be a safe global choice for all simulations, 
we choose p = 1.0002 for our CP® simulations, giving an acceptance rate of 
~ 62%. 


For A updates (where k = 1) the behaviour of p is qualitively very dif¬ 
ferent. For \F^\ > 1.122 we can analytically determine the best value for 
p] to the left of the cusp in Fig. 4.1 (a) the maximum of Eq. (33) occurs 
at the 9 = 7T boundary. To the right the maximum cannot be analytically 
determined and hence we simply settled for a flat value of = 1.28, although 
this choice is of course slightly sub-optimal. Using this method we obtain 
an acceptance rate of ~ 67%. We also tested the possibility of using the 
’’optimized cosh” method for generating random U{1) numbers for the A up¬ 
dates [IS]. We found that in general this method had an acceptance rate of 
~ 90%; however because of the increased numerical cost of the latter method, 
the former gave us a marginally better overall performance. 


4-2. Numerical Errors 

Another issue with the over-heat bath algorithm arises occasionally in the 
calculation of 9oid from the angle between our 0 and T0 vectors. From the 
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scalar product of the two vectors we can calculate cos 9oid to high precision; 
however when cos 9oid is close to 1 , the resultant angle 9oid and ultimately 
sin 9oid cannot be determined to as many signihcant digits because of floating 
point rounding errors. Therefore extreme values of 9oid may be affected by 
these numerical errors. Proceeding with the calculation of (pnew^ the second 
term in Eq. (26) is then likely to introduce a (potentially signihcant) numer¬ 
ical error. In some rare cases 9oid is evaluated to be zero, which would lead 
to a hard failure were the algorithm not modihed. We therefore implement 
the following changes to ensure that the algorithm satishes detailed balance 
to a higher level of numerical accuracy. 

Our solution to the problem is to ’realign’ the vector 


- COS9oid-^'] . (35) 

V \F^\J sm9oid 

when we measure 9oid < 0.001 by subtracting off the numerical error. To do 
this, we take the scalar product of this vector with to determine the 

numerical deviation, i.e.: 


\ = I0-L| cos 9err- (36) 

We then determine a new perpendicular to by hnding: 

- I0±| COS 9err (37) 

Finally we must renormalise this vector to have unit norm. Through our 
simulations we monitored that this method does indeed produce a vector 
perpendicular to F^ up to double precision. 

We still have an issue with the above method when we measure cos 9oid 
to be 1 within machine precision. In this case we determine sin 6 ^ 0 ;^ = 0, 
and so our attempt at constructing a vector perpendicular to F^ would fail. 
The over-relaxation term vanishes in the limit 9oid —t 0; consequently the 
update will either have no effect or change the sign of (j) (on the condition 
we normalise (j) to maintain = 1). There is zero probability of the update 
moving the vector out of this (anti-)aligned state and thus detailed balance 
is broken. We solve this issue by generating a new random unit vector that 
is perpendicular to F^ as our (j)j_ term when choosing the remaining degrees 
of freedom in 0 . 
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5. Finite Volume Effects 


In order to study CSD we measure the integrated autocorrelation times 
of our observables on a series of lattices scaling towards the continuum (by 
increasing /3), while holding the physical volume constant (keeping L/^g 
hxed). Existing studies into the hnite volume effects in the CP® model show 
that L/.^g ^ 10 is a good choice for obtaining sub-percent hnite volume 
effects [H]; however given the high statistical accuracy necessary for our 
simulations we must necessarily extend these studies to understand the hnite 
volume ehects to a greater level of precision. 


In Figs. 5.1, 5.2 and 


5.3| we present the results of our hnite volume 
m, ic and Xt respectively. We 


analysis for the CP® model with f3 = 0.8 for y, 
simulated various values of L, using 50 million conhgurations for the smallest 
values of L, ranging up to 140 million conhgurations for the largest values 
in order to obtain good resolution of hnite volume scaling behaviour. For 
our studies we would like to choose our parameters such that hnite volume 
ehects are under control. Each of our hnite volume scaling plots appears to 
converge to an asymptote within statistical errors (as one would expect from 
the 1/N expansion [H]); it appears that hnite volume ehects for ^g become 
unnoticable only for L/^g ^ 25 with our level of statistics. Interestingly, 
the scaling of ^g appears to show a turning point around ^g — 12 that is 
not predicted by the hrst two terms of the large-iV expansion. Finite volume 
ehects for all other observables appear to be negligible within statistical errors 
from L/^g — 15. 

While the hnite 


In Fig. 5.4 we show the hnite volume scaling of 


volume ehects on this dehnition of the correlation length appear to be more 
under control, the statistical error on the quantity is very large in compar¬ 
ison to those on ^g- H is more accurate therefore to simply use ^g given 
that we hold the ratio L/E,g approximately constant. Where appropriate 
we can use an interpolating function to compensate for mistunings in the 
physical volume L/^g- We remark that the hnite volume scaling behaviour 
is independent of UV ehects and thus these results hold for higher values of 
(3. Lastly we note that ^g does not reproduce the inverse mass gap of the 
theory in the continuum limit. However in the scaling region it should scale 
proportionally to and thus it is valid to use ^g as the dehnition of corre¬ 
lation length when we test the scaling of integrated autocorrelation time for 
our CSD studies. 

Finally we consider hnite volume ehects in the determination of r itself. 
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Figure 5.1: Finite volume scaling for Xm in the CP® model at /3 = 0.8. 
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Figure 5.2: Finite volume scaling for in the CP® model at /3 = 0.8. 
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Figure 5.3: Finite volume scaling for xt in the CP® model at /3 = 0.8. 



Figure 5.4: Finite volume scaling for in the CP® model a.t (3 = 0.8. 
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Figure 5.5: Finite volume scaling for (a) and (6) in the CP® model at /3 = 0.8. 
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Figure 6.1: Continuum scaling of the dimensionless quantity Xt^c- 


In Fig. 5.5 we present the finite volume scaling of (a) and (b) The 
integrated autocorrelation time itself should vanish in the limit L/^g ^ oo; 
however again we use an interpolating function to compensate for any slight 
mistuning of L/^q. 


6. Critical Slowing Down 


In order to study CSD we performed extensive simulations over a wide 
range of /3, with L tuned to give L/^g — 15. This choice ensures that to 
our level of statistical precision we can essentially ignore most finite volume 
effects in the measurements of our observables caused by small fiuctations 
in the quantity L/^g- The numerical results of our simulations are given in 
Table Ib'si 

Before we discuss our results of CSD, we can first use our measurements 
to make a continuum extrapolation of the quantity Xt^c- The action Eq. ([^ 
we employ is correct up to 0{a^) cutoff effects, which is corroborated by the 


apparent linear trend in our scaling plot of Xt^c Tig. 6.1 We can use the 
plot to make a continuum extrapolation of this quantity, to obtain 


XieG = 0.01719(10)(3). (38) 

The central value and first (statistical) error is obtained from fitting the 
correction on ^g > 6.6. The second is a systematic error, which we quote 


15 









/S 


—rem 

Km 

^siow 

Km 

^Xt 

0.8 

6.95(2) 

- 

- 

28.44(9) 

0.85 

10.25(5) 

- 

- 

107.6(7) 

0.9 

15.26(12) 

- 

- 

452(6) 

0.95 

26.86(58) 

20.58(34) 

6.20(54) 

2338(60) 

0.96 

30.48(62) 

22.58(33) 

7.37(56) 

3173(86) 

0.97 

35.30(89) 

24.74(14) 

10.5(1.3) 

4490(150) 

0.98 

39.1(1.2) 

26.54(30) 

11.75(91) 

6120(230) 

0.99 

48.0(1.6) 

29.18(44) 

18.9(1.1) 

9230(380) 

1.0 

61.8(2.1) 

31.96(34) 

30.5(1.9) 

11980(550) 

1.01 

71.1(2.7) 

35.27(22) 

36.5(4.1) 

17220(850) 

1.02 

97.8(3.4) 

38.02(25) 

62.7(4.7) 

27300(1400) 

1.03 

120.1(4.8) 

41.80(22) 

77.5(6.3) 

36800(1800) 

1.04 

180.6(7.4) 

45.54(26) 

139.6(9.0) 

56500(3200) 

1.05 

203.8(9.2) 

49.88(34) 

158(12) 

70900(4500) 


Table 6.2: Summary of integrated autocorrelation times for our measured observables. For 
(3 < 0.9 the slow mode contribution to Xm could not be isolated. 


L 

P 

Stat 

E 

ic 


Xm 


72 

0.8 

80M 

0.6670232(7) 

4.5992(12) 

4.718(20) 

28.0595(18) 

97.03(11) 

96 

0.85 

80M 

0.6222715(5) 

6.3926(20) 

6.60(3) 

46.863(4) 

46.24(11) 

136 

0.9 

80M 

0.5838365(3) 

8.815(4) 

9.07(5) 

78.202(8) 

23.17(11) 

184 

0.95 

lOOM 

0.55026689(20) 

12.095(6) 

12.40(6) 

130.707(15) 

11.96(11) 

192 

0.96 

120M 

0.54404507(17) 

12.869(6) 

13.15(4) 

144.880(16) 

10.66(11) 

208 

0.97 

120M 

0.53797267(16) 

13.709(7) 

14.07(5) 

160.601(20) 

9.53(12) 

224 

0.98 

120M 

0.53204234(15) 

14.597(8) 

15.01(7) 

178.13(2) 

8.05(12) 

232 

0.99 

120M 

0.52624990(15) 

15.526(9) 

15.97(8) 

197.48(3) 

7.12(13) 

248 

1.0 

160M 

0.52058951(13) 

16.528(9) 

16.87(7) 

219.02(3) 

6.48(11) 

264 

1.01 

200M 

0.51505639(11) 

17.593(10) 

18.15(10) 

242.96(3) 

5.67(10) 

288 

1.02 

300M 

0.50964608(9) 

18.721(10) 

19.28(7) 

269.50(4) 

5.09(10) 

304 

1.03 

450M 

0.50435393(8) 

19.934(10) 

20.46(6) 

299.10(4) 

4.34(8) 

320 

1.04 

500M 

0.49917929(8) 

21.206(12) 

21.71(6) 

331.83(4) 

4.01(8) 

344 

1.05 

500M 

0.49410872(8) 

22.558(13) 

23.19(9) 

368.33(5) 

3.40(9) 


Table 6.3: Results of our Monte Carlo simulations of the CP® model. Errors were computed 
using a jackknife analysis. 
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Figure 6.4: Plots of the scaling of (green, below) and (red, above) with ^q. 


as the change in the continnnm value when we ht our full data set while 
accounting also for an term. The value we obtain is in agreement with 
previously measured values for this quantity (e.g. see M), which exhibits 
a signihcant difference from the large-expansion result |15j . 


Xtic — 


1 

2'kN 


0.06 


+ o 



7V=10 


0.153. 


(39) 


To quantify the CSD of our observables we ht the scaling of the integrated 
autocorrelation time with to the expected relations [5], i.e. power law 
scaling for quasi-Gaussian modes, 


r 


ae^ 


(40) 


and exponential scaling for topological modes, 

r ~ a exp {h^^). 


(41) 


In general we expect 6 ~ 2 for quasi-Gaussian modes, although using our 
over-relaxation algorithm we may observe 1 < 6 < 2 [T6] . 


In Fig. |6.4| we display the results of our autocorrelation analysis, where 
we have used an automatic windowing procedure to determine the integrated 
autocorrelation times. The exponential relation for the topological modes 
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Figure 6.5: Plot showing the correlation coefficient r between xt and Xm for the highest 
values of measured. Errors were computed using a jackknife analysis. 


fits very well in the region > 12.09, with y^/dof ~ 1.03 and with htted 
parameters a = 0.1(1.6) x 10“^, b = 12(10), c = 0.24(13). For comparison, 
power law scaling for topological modes is also a reasonable £t, albeit with 
a higher value of x^/dof ~ 1.31, with fitted parameters a = 2.2(4) x 10“^, 
b = 5.56(6). The p-values for these fits are p = 0.41 and p = 0.23 respectively, 
and thus either model gives a valid description of the data. 

However there is a very clear deviation from the expected power law scal¬ 
ing for Xm- We suspect that this is due to a small coupling to the topological 
modes. This idea was suggested in Ref. |6], although the deviations from 
power scaling observed in this study could not be distinguished from mo¬ 
mentum cutoff effects. Here however we take care to £t starting from much 
higher values of furthermore we explicitly verified that also including 


leading-order cutoff effects, i.e. powers of 1/^g ^9^- (^0) and (41), did 

not signihcantly change our results. As a simple test of the coupling between 
the two observables, we compute the correlation coefficient between Xt and 
Xm using 

Cov{Xt,Xm) 


riXuXm) = 


(42) 


A/Var (xt) Var (xm) 

We see very clearly that 
there is a small but highly signihcant correlation between the two. Ultimately 
we are left to consider an alternative ansatz for the scaling of the magnetic 


The results of this analysis is shown in Fig. 6.5 
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Figure 6.6: Plots of the autocorrelation function for ^ = 1.0 data, sampled every 50*^ 
timeslice, before and after the subtraction of the slow modes. The exponential decay of 
the slow modes is indicated by the best ht line. 


susceptibility. 

As originally noted in Ref. [6], the autocorrelation function for the mag¬ 
netic susceptibility displays a long, slow decaying exponential tail that con¬ 
tributes signihcantly to the integrated autocorrelation time at high values of 
the correlation length. We therefore separated the integrated autocorrela¬ 
tion time into two parts. First we take the contribution of the slow modes, 
which we estimate by htting the tail of the autocorrelation function to a 
single exponential decay in a region t ^ t. We were careful to £t the tail in 
a region such that the £t is stable with respect to the time at which the £t 
begins. We then subtract this single exponential mode from the entire auto¬ 
correlation function and then recompute using the standard windowing 
procedure. In the Appendix we provide a brief proof of the validity of this 
method. Our expectation is that the slow modes should exhibit the same 
scaling behaviour as the toplogical modes. We display plots of the scaling 
trends for the slow modes and the remaining modes in Fig. |6.7[ Fitting 


on > 12.8, we can immediately make the following observations: hrstly 
the scaling of the integrated autocorrelation time of the remaining modes 
satisfies the expected power law scaling, with y^/dof ~ 0.58 and p = 0.80. 
We obtain the fitted parameters a = 0.61(2) and b = 1.41(1); the exponent 
b is consistent with quasi-Gaussian scaling [3]. Secondly, the contribution of 
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Figure 6.7: Plots of the scaling of with for both the slow modes and the remaining 
modes after the slow modes have been isolated. 


the slow modes overtakes the contribution of the remaining modes at around 
— 17.1. Lastly, if we £t the slow modes to a power law £t we obtain an 
acceptable £t to the data with x^/dof ~ 1.66, p = 0.11, a = 3.6(1.6) x 10“® 
and b = 5.66(15). Crucially we notice that the exponent b of the power law 
scaling of this mode is in strong agreement with the power law scaling expo¬ 
nent of the topological modes (6 = 5.56(6)); the low value for a is consistent 
with the claim that the coupling to the topological modes is rather weak. An 
exponential £t for the slow modes also gives a £t that is narrowly acceptable, 
with x^/dof ~ 1.93, p = 0.07, although the errors on the htted parameters 
are too large to draw any meaningful comparisons. Finally as a consitency 
check we confirmed that adding the integrated autocorrelation times of the 
slow modes and the remaining modes reproduce results that are compatible 
with our originally measured values. Altogether this analysis suggests that 
a reasonable ansatz for the scaling of the magnetic susceptibility is a double 
power law £t of the form 

r ~ (43) 


(or equivalently the sum of a power law and exponential fit). 

To corroborate our results we also ran a set of simulations with L/^g — 10. 
In Fig. |6.8| we display a plot of our results for the scaling of We remark 


that the scaling of the integrated autocorrelation time again deviates from the 
expected power law trend, and is even more severe in this smaller volume. We 


20 










Figure 6.8: Plots of the scaling of with for our results with L/^a — 10 and 
L/^C — 15. The scaling is fitted to a double power law ansatz. 
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Figure 6.9: Finite volume scaling of the correlation coefficient between xt and Xm in the 
CP® model at /? = 0.8. 
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attribute this to the fact that the coupling to the topological sector appears 
to be a fairly strong hnite volume effect, and thus varies signihcantly as L/^g 
changes. We remark that the hnite volume dependence of was measured 
at a low value of (3 where quasi-Gaussian scaling is dominant. At this value of 
[3 the slow modes that couple to the topological sector can not be isolated and 
extracted as explained above and thus we do not have sufficient information 
to be able to account for hnite volume ehects in the scaling of the slow modes. 
This may acconnt for the observed small deviations from the expected trend, 
particularly for large values of ^g where the slow modes provide the majority 
of the contribution to However we can still measure the correlation 


between Xt and Xm for this data set; a plot of which is shown in Fig 6.9 


This plot demonstrates that the scaling of the correlation coefficient towards 
the inhnite volnme limits appears to be dominantly power-like. Importantly 
we note that for the volumes we can feasibly simulate the correlation is 
non-negligible and thus could introduce a systematic bias on Xm were the 
topological sectors not adeqnately sampled. 


7. Conclusions 

Through our high-statistics simulations of the CP® model we can draw 
several conclusions. Firstly, we have determined that for the level of precision 
reached in our simulations, taking L/^g ^ 15 in the measurement of the 
topological and magnetic susceptibilities is sufficient in order to be able to 
neglect hnite volume effects. However for ^g we can resolve the hnite volume 
ehects up to — 25. 

Secondly we have shown that the CSD of the magnetic susceptibility is 
exacerbated by a small coupling to the topological sector. This results in a 
deviation from the expected power law form of CSD towards the character of 
the CSD of topological modes. Our simulations indicate that this dependence 
is a hnite volume ehect, as the correlation between topological observables 
and other observables (and by extension the CSD of these other observables) 
decreases with increasing physical volume. Nevertheless it follows that it 
is necessary to properly sample the topological sectors to avoid introducing 
systematic biases in the measurement of certain observables. Importantly 
there is no reason for such a feature to be absent in QCD simulations; in 
particular this may have signihcant implications for simulations which are 
frozen to a single topological sector. 
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Appendix A. Analysis of autocorrelation functions 


In this section we present a theoretical description behind the subtraction 
of the slow modes of the autocorrelation function for the magnetic suscepti¬ 
bility. We assume that Xm decomposes into two parts, 

Xm = Xm + CXt, (A.l) 


which follows from the correlation we measured between these obserables. 
The quantity x'm represents the contribution of modes which are not corre¬ 
lated with Xt- 

We insert this expression into the formula for the autocorrelation function 
Eq. (16), thus obtaining 


+ 2Cov^^^^j(t), 


(A.2) 


where represents the cross-correlation between the observables. 

We neglect any remaining cross-correlation between these two observables, 
i.e. we take = 0. 

Lastly we assume that the autocorrelation functions display exponential 
asymptotic behaviour: 

~ ai exp (A.3) 

V WL / 

~ “ 2 exp ■ (A-4) 
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It follows therefore that when we have S> & 2 /'rxt fhe exponential be¬ 

haviour of C^^{t) will provide the dominant asymptotic behaviour of 
At large enough values of (3 this inequality will be satished; htting the tail 
of the autocorrelation function at asymptotic times will thus yield the con¬ 
tribution of the slow modes to the integrated autocorrelation time. 
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